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I. INTRODUCTION 


Accurate path control of surface ships and underwater vehicles along prescribed 
geographical paths is a basic problem that is becoming increasingly important, par- 
ticularly as the missions of ocean vehicles become more complicated with strict 
requirements for performance. In order for a control law to be able to perform its 
mission in a realistic operational scenario it has to be robust enough so that it can 
maintain stability and accuracy of operations in the presence of modeling errors and 
environmental uncertainties. The robustness properties of the design are particu- 
larly important due to the unpredictable nature of the ocean environment and the 
changes in the hydrodynamic characteristics of the vehicle during turning, changes 
in the forward speed, or operations in proximity to other objects in the area. For 
these reasons, there exists a need for the analysis of the robustness characteristics 
of a particular control law design and the establishment of a rational operational 
envelope based on stability and performance criteria. Previous studies [Parsons 
and Cuong (1977)] showed that gain adaption is highly desirable due to changes 
in the linearized vehicle hydrodynamics with different operation conditions, such 
as depth under keel. The resulting adaptation scheme [Parsons and Cuong (1980)] 
required significant vehicle motion, which may be undesirable when operating in 
restricted waters, or in object recognition and localization tasks. Integral control 
techniques [Parsons and Cuong (1981)] proved quite effective, but neglected the 
nonlinear behavior of the vehicle, which becomes very important at low speeds and 
hover operations. Model based compensators exhibit robust behavior under con- 


ditions of parameter uncertainty, which is as good as the classical linear quadratic 


regulators for linear output feedback systems [Healey (1992)]. Alternatively, sliding 
mode controllers exhibit very robust characteristics given an estimate of the param- 
eter uncertainty and/or disturbances [Papoulias and Healey (1992)], [Yoerger and 
Slotine (1985)]. Sliding mode control, however, does not offer an infinitely robust 
design and it suffers from a series of bifurcation phenomena and loss of stability 
unless proper care is exercised [Papoulias (1991)]. 

In this work we analyze the problem of loss of stability of a marine vehicle 
under cross track error control in the presence of mathematical versus actual system 
mismatch. For demonstration purposes, variations in the heading angle control 
gain are studied. Previous studies [Oral (1993)] concentrated on system response 
assuming perfect and complete state measurement. Particular emphasis in this work 
is placed on analyzing the effects of observer design on system response after initial 
loss of stability of straight line motion. The main loss of stability cases analyzed 
here occur in the form of generic bifurcations to periodic solutions (Guckenheimer 
and Holmes (1983)]. We use center manifold reduction techniques and averaging 
in order to capture the stability properties of the resulting limit cycles [Chow and 
Mallet-Paret (1977)]. It is shown that the dynamics of the observer may have a 
significant effect on the computed gain margin of the control system depending 
on the particular basis used. All computations in this work are conducted for the 
NPS autonomous underwater vehicle [Bahrke (1992)] and all results are presented in 
standard dimensionless quantities with respect to vehicle length, 7.3 ft, and nominal 


forward speed, 2 ft/sec. 


II. PROBLEM FORMULATION 


A. EQUATIONS OF MOTION 
The linear maneuvering equations of motion of a marine vehicle in the hori- 


zontal plane are written in dimensionless form as, 
miv+r+z26r)=Y;r+Y,v+Y,r+Yuvt+ Ysé (27) 


Lr+mzg(vt+r)=Nyrt+ N,v + Nr + Nv + N56 (a) 


where all symbols are explained in the nomenclature. Equations (2.1) and (2.2) can 
be used to derive a second order transfer function between the rudder angle 6 and 
yaw rate r. For low frequency maneuvering motions this second order equation can 
be approximated by expanding in Taylor series and keeping the first order terms 
only. The result is 

r=ar+ 06 (2.3) 


Equation (2.3), which is sometimes referred to as Nomoto’s first order model, is 
particularly useful in control system design since no sway velocity feedback is neces- 
sary. This equation predicts linear variation of the steady state turning rate versus 
rudder angle. In reality, the r-6 curve has characteristics of softening spring mainly 
due to speed loss during turning. To account for this a modified version of equation 
(2.3) is used, 

r= ar+agr° + 66 (2.4) 


where a3 is usually determined from steady state results. Finally, the model is 


complete by the incorporation of the kinematic equations, 


W=r (225) 


y=sin (2.6) 


where W is the vehicle heading, and y is the cross track error off a desired straight 


line path. 


B. COMPENSATOR DESIGN 

In control theory it is known that the eigenvalues of the controller are not 
affected by the eigenvalues of the observer. This allows us to design the controller 
and observer separately which is known as the separation principle. The combination 
is called a compensator. 

Equations (2.3), (2.5), and (2.6) govern the steering control of the model used 


in this section. The control law can be expressed as, 


5 = dsa tanh (=) (25) 
beat 


where around the nominal state VY = r = y = 0 we have 
bdo = KyV + K,r+ Kyy (2.8) 


6 is the rudder angle and Ky, K,, and Ky, are the control gains of the system. 
The linear control law is d9. The rudder angle 6 is defined by a hyperbolic tangent 
function to include the saturation to our problem as shown in Figure 2.1. Saturation 
occurs at és, which is the saturation limit generally taken as 0.4 rad. 


The linearized form of equations of motions in the vicinity of V =r = y = 0 


are, 
W=r (2.9) 
r=ar + bédo (2.10) 
y= (2.11) 





Figure 2.1: Saturation in 6. 
These equations can be expressed in state space form as 


X=AX+4+ Bu (212) 


where 


is the state vector, 


is the control distribution vector. 


The observer equations are 
Pe = AVG esa L(Y CX) 213) 


where X is the estimate of X, Y is the output of the system Y = y, and C is the 


sensor vector C' = [0 0 1}. 


co A | 


The error in the estimate of X is defined by 
Xa (2.14) 
Using equations (2.12), (2.13), and (2.14) we can obtain 


X =(A—LC)X (2.15) 


We can rewrite equations (2.9), (2.10), and (2.11) in the form of 

Ny O10) fu 0 

; Oa od] tri+lols (2.16) 
y Pn) y 0 


Y 
Y =(00 1] | r | (2.17) 


y 


ly 
ie | e, | (2.18) 
ty 


After performing the matrix operations we obtain 


0 1 —by 
=|0a —é, 
10 —-@, 


Using equation (2.13) we can rewrite equation (2.8) as follows, 


X= 





and 


The observer gains are, 


oe 
as 


X= 


Cae “ZJre Heit 
2 ie fa 1 


| 219 


$9 = Ky(W — ¥) 4+ K,(r -7) + K,(y -9) 


Finally, we can write our compensator equations in the form 


Y 010 0 0 0 N 
r 0 a 0 beKy 6K, bK, r 
Sie eo of oo y 
xy! I|oul 1000 O I le v 
000 0 a -£ ; 
bah ho 00 1 08 HORE 











If we look at the matrix carefully we will see that it is in the form 
X|_[A-BK BK x 
je) = 0 A-—LC X 
which has the following characteristic equation, 
det(A — BK — sI|det{A -—LC — sl] =0 


This indicates that the dynamics of the observer are completely independent of the 


dynamics (eigenvalues) of the controller. Thus A and L can be designed separately. 


C. CALCULATION OF CONTROL GAINS 


A is the Jacobian matrix of the system 


0 1 0 
Az=| bKy a+bK, bK, 
1 0 0 








The characteristic equation of the matrix A is 
* — (a + bK,)r? — bK yA — OK, = 0 
If the desired characteristic equation has the general form 
A + agr*? + aA + a0 = 0 


the control gains can be found as 


rahe | 
oa ; 
= Qgat+a 
—— ; 
a 
K, = ary 


The desired characteristic equation can be written with respect to the desired natural 
frequency and some optimum coefficients. The ITAE criterion for a third order 
equation is 


3° + 1.75w,s? + 2.15w2s + w® = 0 
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where w,, 1s the desired controller natural frequency. 


Therefore the contro] gains can be calculated for a given natural frequency, as 


a = 2.15w? 
Q2 = 1./5w, 
ao = w> 


D. CALCULATION OF OBSERVER GAINS 


If we define A as the Jacobian matrix of the system 


01 —by 
A=|0a —4, 
10 -2, 


the characteristic equation of the matrix A is 
\* + (£, — a) + (€y — aly)rA + (€, — aly) = 0 
If the desired characteristic equation has the general form 
A? + 2A? + WA +o = 0 


the observer gains can be found as 


Ey = a oa 2 
ly = al, aa 1 
fl, = aly+o 


Applying the ITAE criteria, observer gains can be calculated for a given natural 


frequency as 
71 = Dp 15w?, 
Yo = 1.75wWnpo 


Yo = Who 
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where Wag 1s the observer natural frequency. 


E. CHARACTERISTICS OF ITAE CRITERIA 

In the calculations of gains we applied the ITAE criteria. If we look at Figure 
2.2 for the step response of ITAE, we see that the response gets faster as the natural 
frequency increases. For example, the settling time is 10 normalized seconds, or 10 


seconds for w, = 1, 1 second for w, = 10, and so on. 


STEP ResrOnNsiesot 1A 
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Nommalized Time ( t*wn ) 


Figure 2.2: Step response of ITAE. 
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Ill. HOPF BIFURCATION 


A. INTRODUCTION 

An important quantity in assessing the robustness of a particular control law 
design to parameter variations and unmodeled dynamics is the gain margin. This is 
defined as the extent to which changes can be inflicted on the system gain without 
loss in stability. to this end, we assume that the heading error gain Ky is multiplied 
by a constant C’. By definition, of Hopf bifurcation occurs when a pair of complex 
conjugate eigenvalues cross into the right hand half-plane. When this occurs the 
system will deviate from a steady solution in an oscillatory manner. This deviation 
can be either supercritical or subcritical [Seydel (1988)]. As the parameter C' crosses 
the critical value, one pair of complex conjugate eigenvalues of the linear system 
matrix crosses transversely the imaginary axis. Locally, as C approaches C',;,, the 
periodic solutions are located on the two dimensional Eucledian plane spanned by the 
eigenvectors of the Jacobian matrix of the system which corresponds to the critical 
pair of eigenvalues. In this chap ~- stability properties of the periodic solutions are 
established. In order to establish those properties the main nonlinear terms that 
dominate the system are isolated. Center manifold theory is used to reduce the 
flow to a two dimensional manifold. The method of averaging is then applied to the 
reduced system. 

The critical value of c for stability of straight line motion remains the same as 
[Oral (1993)], which is 

Cerit = 0.2658 
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This is because the dynamics of the controller are independent from the dynamics 


of the observer as explained in Chapter II. 


B. THIRD ORDER EXPANSIONS OF THE SYSTEM EQUATIONS 
1. Perturbation in Ky 
In the previous chapter we worked on the linear system. Now we are go- 
ing to introduce the nonlinear terms to our compensator. In this case the equations 


of motion are 


wo =r (3.1) 
r = ar+agr° +66 (3.2) 
y = sin (3.3) 
where 
bo 
6 = dca: tanh (=) (3.4) 
beat 
5 = CKy(U—WV)4K,(r—F)+K,(y —9) (3.5) 
or in compact form, 
“ ig 
XSF (2) — | ee a (3.6) 


This system can be written in the form 
X = AX +4(z) (3.7) 


A is the Jacobian matrix of f(z) evaluated at X = 0, and g(z) contains all nonlinear 
terms of Equations (3.1), (3.2), and (3.3). Taylor expansion of the nonlinear terms 


about the equilibrium, where we keep the first non-vanishing nonlinear coefficients 
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only, gives 


snV = W- -¥° (3.8) 
1 3 
6 = do a 362, 00 (3.9) 


Substitution of Equations (3.8) and (3.9) into Equations (3.1), (3.2), and (3.3) gives 


us the A matrix in Equation (3.7) as follows, 





0 1 0 0 0 0 
bcKy a+bK, 6K, —bcKy -—bK, —bK, 
1 0 0 0 0 0 
aS | ip 0 oF oO ee (3.10) 
0 0 0 0 a —. 
0 0 0 1 te 
The nonlinear parts are, 
0 
b 
a3r : 362, 5: 
g(z) = ae (3.11) 
0 
0 
0 


If we introduce the transformation matrix (T) of eigenvectors of A evaluated at the 


bifurcation point, 


T [m,,] t,J = 1,2,3,4,5,6 (3.12) 


/ aie [n;;] peels 253545586 (sal) 


the linear change of coordinates, 
coal ae z=T's (3.14) 
transforms system (3.7) into its normal form 


= TATz + T7'9(Tz) (3.15) 
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With this transformation we get, 


0 —Wo 0 0 0 0 

w 0 0 0 0 O 

|| hee OP OOO 
a 0 0 O PF 0 0 
0 O O 0 Ps 0 

tO O 0 0 0 F 


(3.16) 


Using center manifold theory we get, as shown in [Chow and Mallet-Paret (1977)], 


0 
bo, 22 + boo z2 + bogzeze + boazr rf 
3127 + l30z5 + bg3ziz2 + bs4z125 
0 
0 
0 


[el 


Substitution of Equation (3.14) into Equation (3.7) yields, 


3 2 2 3 
Zp = ~—Wo22 +7112) +7122] 22 + 7132122 + 71422 


Z2 


2. Integral Averaging 


We write Equations (3.18) and (3.19) in the form 


21 —wo2z2 + Fi(2, z2), 


Z2 W021 + Fo(21, 22) 
If we introduce polar coordinates in the form, 
z, = Rcos8, za = Rsin#d 


Equations (3.20), (3.21) result in 


R 


F,(R, @) cos @ + F(R, @) sin @ 


RO woR + Fo(R,9) cos 6 — F(R, 9) sin 8 
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3 2 2 3 
W021 + 7212] + 7222, 22 + 1232122 + 172429 


(3.17) 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


(3.22) 


(3.23) 


(3.24) 


Equation (3.23) yields 
R= P(6)R° 


(3.25) 


where P(6) is a 27-periodic function in the angular coordinate 0. If Equation (3.25) 


is averaged over one cycle in 8, we get an equation with constant coefficients, 


R= KR 
where, 


] nr 
K = = | P(6)- do 


Equation (3.27) is simplified after evaluation of the integral as, 
1 
K = 8 [3ri1 + 713 + 722 + 3724] 


where the coefficients are as follows, 


Tyy = Niza + rises 
M13 = Ni2log + 213034 
[22 = No2lo3 + No3l33 
[24 = Nogloo + No3l32 


(3.26) 


(3.27) 


(3.28) 


where the nj2, 213, N22, and n23 are the elements of the inverse of transformation 


matrix 7’. The values of the coefficients @;; are given by the following expressions 


toy = apm, — by | Kgms, + Kems, + Kom}, + 30° KG K-miyms1 


+3c? Ki. Kym?,me) + 3cKy K?m2ima Te 3cKy Kimsmai + 3K, Kims.msi 


43K? Kym? me, + 6cKy K, Kymamsimey| 
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lo3 


a3m3, — by lo Kgms, + Kims, + Kimé, + 3c? Ky K,mi.ms2 

+3c* Ki, K,mi.me2 + 3cKy K?mi,m42 + 3cKy Kims.ma2 + 3K, Kimé.ms2 
+3K? Kymj.me2 + 6cKy K, Kymams2meo| 

3a3m3,™M22 — by [Sc K gm, ma + 3K3m?2,ms2 <5 3Kimé, M62 

$30? KK, (miymse + 2maimaams) + 3c°KG Ky (m3, mez + 2m4ma2me1) 
+3cKy K? (m2,ma2 + 2msims2M41 ) ~ 3cK y Ks (ma,ma2 + 2meime2ma1) 
+3K, Ke (mamse 5 2me1m62™s1 ) + 3K°K, (m3,me2 + 2msims2me; | 
+6cKy K,K, (71741™M51M¢62 ate (1141™M52 1 mM42M51) me1)| 

3a3m3,m3. — by [3c°.K 3 mam, +3K?msim3, + 3K}meimé, 

+3c°K ai, (m3.ms5, + 2maM42™52) + 3c°K he (mi.me a 2mam42mey)] 
+3cKy kK? (m2oma1 oa 2ms1M52M42 + 3cKy ne (maomat == 2mermezma2) 
+3K,K? (mg.msi + 2meime2ms2) + 3K2K, (m3mer + 2msimsomen) 


+6cKyK,K, (m42ms2me1 + (™41™M52 + M42M51) me2)| 


l 
fz, = —emi, 
l 
f32 = emis 
I 2 
C35 = 9M 12 
I 2 
ls, => ~ 5 M11 12 
b 
a 362, 
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C. RESULTS 

The value of AK is important for us to determine whether the bifurcation 1s 
supercritical (A < 0) or subcritical (K > 0). In this study we wanted to see the 
effects of observer dynamics to our system, especially the value of K. To do that 
we used the Fortran code (Appendix A) for the numerical results. The result was 
that the value of K was not affected by changes in the observer natural frequency. 
The reason for this can be traced back to the definition of K, Equation (3.28). It 
can be seen that in the definitions for r;; and @;; only the first two eigenvectors mj, 
m2 fori = 1,2,...6 of the A matrix, Equation (3.10) appear. Therefore, we have 
to show that these remain the same for all observer natural frequencies. 

This A matrix, Equation (3.10), actually consists of 4 block matrices, each 


3 X 3, which are the same as shown in Equation (2.22). Let us denote the A matrix 


as follows, 
A B 
a=[9 | 
The eigenvalues of A can be computed by 
A— XI B ~9 
0 C—rAI | 


or 


[A =| eT | 0 


We can group the eigenvalues in two different groups: \,; for 7 = 1,2,3 are the 
eigenvalues of A an A2; for : = 1,2,3 are the eigenvalues of C. The eigenvectors 
associated with these eigenvalues can be found as follows. 


For A = ALis 
A — Ayal B Vv) — 0 
0 C anal Ay il v2 7 0 


i 


and 


lI 
= 


[A — Ai iJ] [vi] + [B] [v2] 


ll 
=) 


[O}[oy] + [C — A1,i2] [v2] 


Since 4, 1s an eigenvalue of A and the eigenvalues of A and C are distinct, the 


matrix [C — ,,,J] is nonsingular which means that [v2] = 0. Therefore, we get 
[A = Ai! ][v1] = 0 


which means that v; is an eigenvector of A. Therefore, the first three eigenvectors of 
A are the same as the eigenvectors of A and they are independent of the dynamics 
of the observer. Of course, the remaining three eigenvectors of A depend on the 
oberver natural frequency, but, as we pointed out earlier, none of these appear in 


the definition of the nonlinear stability coefficient K. 
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IV. COMPENSATOR IN A DIFFERENT BASIS 


A. CRITICAL VALUE OF C 

If we look at Equation (3.6), we can see that the basis for our system was 
[X,X]. Now we are going to represent our system in [X,X] basis where X is 
the estimate of X. In this compensator basis the critical value of C in Equation 
(4.3) is no longer constant. Therefore, we used a Fortran code (Appendix B) to 
calculate the critical C’ values for different observer natural frequencies. A plot 
of these critical C’ values for different observer natural frequencies can be seen in 
Fig. 4.1. The observer natural frequencies are in nondimensional seconds whereas 
the control natural frequencies are normalized with respect to the corresponding 
observer natural frequencies. The system is unstable for values of C' less than the 


critical value. 


B. THIRD ORDER EXPANSIONS OF THE SYSTEM EQUATIONS 
1. Perturbation in Ky 

In the previous chapters we worked on the [X,X] basis. Now we are 

going to represent our system in the new basis, which is [X, X ], where X is the 


estimate of X. Our equations of motion were, 


War 
r = ar+ag3r° +66 (4.1) 
y = sin 
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CCRIT vs WN(normalized) 
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Figure 4.1: Cai, versus natural frequency for Ky. 
where 
So 
6 =Widy, tant (4.2) 
on 
6 = CKy¥+K,74+Kyj (4.3) 
or in compact form, 
~ ae: 
PC ile il Ae Wor,y,¥,7, 9] (4.4) 
This system can be written in the form 
X = AX + 9(z) (4.5) 
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A is the Jacobian matrix of f(z) evaluated at X = 0, and g(z) contains all non- 
linear terms of Equation (4.1). Taylor expansion of the nonlinear terms about the 


equilibrium, where we keep the first non-vanishing nonlinear coefficients only, gives 


sn¥Y = ¥--0 (4.6) 
5 oe ae (4.7) 
a3ne ° : 


Substitution of Equations (4.6) and (4.7) into Equation (4.1) gives us the A matrix 


in Equation (4.5) as follows 


01 0 0 0 0 
0a 0 bcKy OK, bK, 
10 0 0 0 0 
aia INOMO. fe c0) 21. See eS) 
00 4 bcKy OK, —é,4+ 6K, 
10 0 4 1 0 —, 
The nonlinear parts are, 
0 
b 
ei 
Sarg 
g(z) = 6 (4.9) 
0 
b 43 
~ 397, °0 
0 


If we introduce the transformation matrix (T) of eigenvectors of A evaluated at the 


bifurcation point, 


T 


[m;,] 2,7 = V2. 3,4,9,6 (4.10) 


T> In:;] i,j =1,2,3,4,5,6 (4.11) 


the linear change of coordinates, 


niawize. ae (4.12) 


nA | 


transforms system (4.5) into its normal form 


z= TAT Fogle) (4.13) 
At the bifurcation point 
0 -w 0 0 0 O 
<1 47 | 2020 Sele oe 0 
TU Tenet (ees Ane (4.14) 
0 0 O O P 0 
“0 0 0 0 0 Pa 


with wo > 0 and P; < 0. 
Using center manifold theory we get, as shown in (Chow and Mallet-Paret 


(1977)], 
0 
boi z7 + b2025 + Logzize + b242127 
3127 + lgo2z3 + lggzpze + b342123 


0 (4.15) 
bs, 2z7 + l5225 + — + 542,23 | 
Substitution of Equation (4.12) into Equation (4.5) yields, 
Z) = —wozetriizy t ri2z22 + 7132124 + M1425 (4.16) 
Zp = Wor taize + reezez2 + 7232125 + 7242p (4.17) 


where the terms r,;; are evaluated by a Fortran code (Appendix C). 
For values of C close to its critical value, Equations (4.16) and (4.17) 


become, 
. as / / 3 2 2 3 4 18 
Zz) = alex — (wot w'e)z2 + riizy + 7122122 + 7132122 + 71422 ~~ (4.18) 
A xa t f 3 2 4 3 
zg = (wot w'e)z + a’ez, +112) + 1222] 22 + 7232125 + 12429 (4.19) 


where ¢€ is the difference between C and its critical value Cy. The terms a’ and w’ 
denote the derivative of the real and imaginary part, respectively, of the critical pair 


of the eigenvalues with respect to C’ evaluated at Cy 


Ze 


2. Integral Averaging 


We write Equations (4.18) and (4.19) in the form, 


Zz) = a’ez, — (wot w'é)z2 + Fi(21, 22) (4.20) 


29 (wo + we) zy +p a’ Ez + F(z, z2) (4.21) 
If we introduce polar coordinates in the form, 


z, = Rcos6, z2 = Rsin# (4.22) 


Equations (4.20), (4.21) result in 


R = a'eR+F(R,0)cos6 + F(R, 6) sin 6 (4.23) 
RO = (wo+w'e)R+ F(R, 6) cos 6 — F,(R, 6) sin 8 (4.24) 

Equation (4.23) yields 
R=a'cR+ P(6)R (4.25) 


where P(@) is a 27-periodic function in the angular coordinate 6. If Equation (4.25) 
is averaged over one cycle in 0 [Chow and Mallet-Paret (1977)], we get an equation 


with constant coefficients, 


R=a'eR+ KR? (4.26) 
where 
l rd 
K == I P(0)d6 (4.27) 


Equation (4.27) is simplified after evaluation of the integral as, 
1 
ie 3 (3ri, + rig + r22 + 3724] (4.28) 


where the coefficients are as follows 
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Ti, = 42021 + N33) + 21525; 


713 = 142l04 + 213634 + Ny5L54 
T22 = Nobo + 22333 + N25e53 
124 = Nolo + N23032 + No5b52 


where the nj2, 213, 215, 222, 223, and n25 are the elements of the inverse of trans- 


formation matrix J’. The values of the coefficients €5;, C22, £23, €24, £31, €32, €33, and 


€34 are the same as in Chapter III. The values of the other £;; coefficients are given 


by the following expressions: 


51 


53 


—by °K yma, + Kem§, + Kemg, + 3° KG K,-miyms 

+3c? Ki, K,m2,me + 3cKyK?mi,m4; + 3cK y Ki mi,ma + 3K, Kimgmsi 
+3K? Kymiime: + 6cKyK,Kymaimsimei| 

—by [PKG mi, + Kem, + Kom, + 3c? KG K,mi,ms2 

+3? KU Kymi.me2 + 3cKy K?m},maz + 8cKy K2miama2 + 3K, K2megms2 


+31? Kymi,me2 + 6cKy K, Kymazms2m60| 


—by [3° KG m3,ma2 + 3K2m2,msq + 3K2m3 mee 

+3c° Ki K, (m2 mse + 2maima2ms; ] + 3° KK, (m2, m62 +. 2maim42mMe; ] 
+3cKy K? (m2, maz + 2msims2™M41 ) + 3cKy Ki (m3, maz + 2merme2Ma: | 
+3K,K? (m2, msz + 2meime2ms1) + 3K? Ky (m2, mez + 2ms1™Ms52™Me1) 


+6cKyK,K, (maimsime2 + (™41™M52 + M42™M51) mei)| 


24 


— 


bs4 = — by 3K y mama, ta 3K?ms1ms, + 3K} mei mé 
+3c°KiK, (mi.ms1 + 2mama2ms2) +3°K2K, (mi,me 4 2m41m42™¢2) 
+3cKyK? (m2.ma - 2msims2mM42) + 8cKyK? (maomar aa 2me1mezMa2) 
+3K,K3 (m3.ms1 + 2merme2Ms2) +3K?K, (m3,mer + 2ms1Ms2™e¢q2) 
+6cKy K,K, (M42™M52™Me, == (™m41M52 = ™M42™s5} ) me2)| 
C. RESULTS 
Existence and stability of limit cycles can be determined by analyzing the 
equilibrium points of the averaged Equation (4.26), which correspond to periodic 


solutions in 2), 22 as can be seen from Equation (4.22). From Equation (4.26) we 


can easily see that: 
1. If a’ > 0, then 


(a) if K > 0, then unstable periodic solutions co-exist with the stable equi- 


librium for e < 0, and 


(b) if K <0, then stable periodic solutions co-exist with the unstable equi- 


librium for € > 0. 
2. If a’ < 0, then 


(a) if K > 0, then unstable periodic solutions co-exist with the stable equi- 


librium for € > 0, and 


(b) if K < 0, then stable periodic solutions co-exist with the unstable equi- 


librium for e < 0. 
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We refer to K < 0 as the supercritical, and AK > 0 as the subcritical PAH bi- 
furcation. In the supercritical case, after the equilibrium state loses its stability 
the system converges to a stable periodic solution with amplitude which increases 
continuously as the difference € is increased. 

In the subcritical case, however, before the equilibrium state loses stability, its 
domain of attraction becomes very small since it is bounded by the amplitudes of the 
unstable limit cycles. In such a case, an initial disturbance of sufficient magnitude 
can throw the system off the nominal path even before its domain of attraction 
has completely shrunk to zero. As the nominal equilibrium becomes unstable, the 
system jumps to a different state of motion with a locally, at ¢ = 0, discontinuous 
increase in the amplitude [Papoulias (1993)]. 

In our case, the value of a’ is always negative, which can be seen easily from 
Figure 4.1. If we look at the nature of the curve for the critical value of C’ for 
different natural frequencies, we will see that as the value of critical C decreases for 
the same natural frequency, the system becomes unstable. 

After using a Fortran code (Appendix C) we observed that the nonlinear sta- 
bility coefficient AK depends on observer dynamics. Figure 4.2 shows that for a given 
control design, the observer must be as responsive as possible to ensure negative K 
(stable limit cycle). On the other hand, for a given observer design, if the control 


law is too slow we get subcritical behavior (unstable limit cycle). 
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Figure 4.2: 


iin ELEC 2 I Oe 
WN(normalized) 


Kx, versus w, for different observer w,. 


bo 
=} 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

An investigation of the nonlinear dynamic response characteristics of a marine 
vehicle has been presented. Particular emphasis in this work was placed on analyzing 
the effects of observer design on system response after initial loss of stability of 
straight line motion. Bifurcation theory techniques were utilized in order to assess 


that behavior. The main conclusions of this work can be summarized as follows. 


1. There exists a critical point for a certain combination of system gains and 
system parameters for stability of straight line motion. The loss of stability 
occurs generically in the form of Poincare-Andronov-Hopf bifurcations. As 
the parameter crosses its critical value, a family of periodic orbits, self sus- 
tained oscillations develops. Center manifold reduction and integral averaging 
techniques were used in order to establish the direction of the bifurcation and 


stability of the resulting periodic solutions [Papoulias, Oral (1993)]. 


2. For [X, X] basis the critical point does not depend on the observer dynamics 
(separation principle). The nonlinear stability coefficient K was not influenced 
by observer dynamics. The previous reduction process shows that A depends 
on the first two eigenvectors of the 6 x 6 matrix A. Matrix algebra shows that 


these eigenvectors are associated only with the controller dynamics. 


3. For [X, X] basis the critical value depends on observer dynamics. For a given 


control design, the observer must be as responsive as possible to maximize the 
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region of stability. On the other hand, for a given observer design, the control 


must be as slow as possible to maximize region of stability. 


4. The nonlinear stability coefficient K depends on observer dynamics for this 
basis. For a given control design, the observer must be as responsive as possible 
to ensure negative K (stable limit cycles). In this benign loss of stability 
the resulting periodic solutions are continuous single-valued functions of the 
parameter distance from its critical value. On the other hand, for a given 
observer design, if the control law is too slow we get subcritical behavior 
(unstable limit cycles). In such a case, the periodic solutions develop with 
what appears to be a discontinuous increase in the amplitude of oscillations 


[Papoulias, Oral (1993)]. 


B. RECOMMENDATIONS 


The differences between the two bases with respect to robustness properties of 


the system have to be analyzed. 
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APPENDIX A 


HOPF BIFURCATION PROGRAM FOR [X, X] BASIS 


23456 


R? >? RK HK HK K| 


PROGRAM HTKPSI 
HOPF BIFURCATIONS 
NOMOTO’S FIRST ORDER MODEL 


CALCULATIONS FOR K AND CCRITICAL IF KPSI CHANGES W/ C 


7 
SIPLICIT DOUBLE PRECISION (A-H,0=Z) 


meebo oe PRECISION K1,K2,K,LPSI,LY,LR,1Z,L, 

MASS, NV,NR, NVDOT, NRDOT, NDRS, NDRB, KPSI, KR, KY, K3, 
ae be, 2e 247d, L353 Sine, Lol, Lo2 bos, L54, 
Mi1,M1l2,Mi3eM14, Mis; Ml om@hM21,M22,M23,M24,M25,M26, 
M31,M32,M33,M34,M35,M36,M41,M42,M43,M44,M45,M46, 
Moi, M52 , Meee Ma4 7 Maer S67, Mow, M62, M63 ,M647M65, M66, 
Pie Nes, NS, Nia NIS> Nic, NownNZ2,N23,NZ49N25,N26, 
N31,N32,N33,N34,N35,N36,N41,N42,N43,N44,N45,N46, 

Nei No, Noo, N>4 Noo, Noo, NGlgiie2 ,N63,N64,N65, N66 


DIMENSION AMAT(6,6),T(6,6) ,TINV(6,6),FV1(6),IV1(6),YYY(6, 6) 
MEMENSTON WR(6) ,WI (6), TSAVE(6,6) , TLUD(6,6) , IVLUD(6) , SVLUD(6) 
DIMENSION ASAVE (6,6) ,A1(6,6) ,A2 (6,6) 

OPEN (11,FILE=’AKPSI.MAT’ , STATUS=’ NEW’ ) 


WEIGHT=435.0 


IZ = 4 500) 

L =e. 5 
RHO =e, 94 

G 8 PR 
XG =r LO4 


MASS =WEIGHT/G 

Meteo §8=60—=MASS/ (0.5*RHO*L**3) 
IZ Mee * RHO*L,* *5) 
XG =XG/L 

memOT =-0.00000 

MYDOT =-0.03430 

YR =) 00000 

YV =-0.10700 

memo =+0.01241 

YDRB <=+0.01241 

NRDOT =-0.00047 

NVDOT =-0.00000 

NR =-0.00390 
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© Ole) @. 


205 


50 


=. 00000 

Org 57 ~ UDns 

NDRB: =+0.28354ORS 

DH =(IZ-NRDOT) * (MASS-YVDOT) - 
(MASS*XG-YRDOT) * (MASS* XG-NVDOT ) 
(IZ-NRBOT) *YV- (M&SS*AG-YRDOT) “NV) 7 BH 
(IZ=NRDOT) *(-MaSS+ YR) — 
(MASS*XG-YRDOT) * (-MASS* XG+NR) ) /DH 
( 

( 

( 

( 


2 


( 
( 


MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
MASS-YVDOT) * (-MASS* XG rhe 

MASS* XG-NVDOT) * (-MASS+YR) ) /DH 

ILZ=-NRDOT) *YDRS= (MASS "XG -vE Pe) NPre yy ba 
BiZ= ( (IZ-NRDOT) *“YDRB=- (MASS? XG-VeReeiDe. | eine) ee 
B21=( (MASS-YVDOT) *NDRS- (MASS* XG-NVDOT) *YDRS) /DH 
B22=( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 
Bisbee 

BZ =e eb 


NoMa 
= ( 


Wie (e704) 

READ Cs ae) WNMIN , WNMAX , IWN 
INCR=IWN 

WRITE (*,*) ‘ENTER OBSERVER WN 
READ (*,*) WNO 

WEE (eer } 

READ WG* joe 

DO is Dsat 


WRITE (*,1006) 

READS (t;”) DO 

WRITE (*,1008) 

READ (>) *) CeraE 
Cl=(A11*A22-A21*A12)* (A214 ei ee 


C2=(A11+A22) * (A21*B1—-Al1“B2)+B2"* (All “A220 -2oa ee 


C3=— (420781 Alle 


BJO) AB abesdl SASS 


WN =WNMIN+ (WNMAX-WNMIN) * (II-1) / (INCR-1) 
jeu e@nbaicme Asi iaal 

ALPHAO0=WN* *3 

ALPHA1=2 .15*WN**2 

ALPHA2=1.75*WN 

KPSI=-ALPHA1/B 

KY =-ALPHAO/B 

KR =- (ALPHA2+A) /B 

GAMA0 =WNO* * 3 

GAMAI=2 -15=ViNOr=2 
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GAMA2=1.75*WNO 
LY=A+GAMA2 

LPSI=A* LY+GAMA1 
LR=A* LPSI+GAMAO 


2234567 


™ 


NUP WNHORHDUP WHR DUP WNHORHDUPWNHR DUP WNF OY UIP W 


r~— ~~ er ir ie ae ess ee ess aes ass ae ee ae aes ae ee ees ae eee ae 


0 
Beemer KPSI 


,6 
,J) =AMAT (I,J) 


eee RG(6,6,AMAT,WR,WI,1,YYY,IV1,FV1, TERR) 
CALL DSOMEG (IEV,WR,WI,OMEGA, CHECK) 
mere (*,*) IEV 


33 


fis 


Zul 


14 


ZZ 


ile 


23 


16 


24 


ey 


WRITE (60,%*) (WR{ ITREAD) fReAn= ie 


OMEGA0 =OMEGA 
DOr > Litre 


T( 1, Daa ee 
T( 1, Zee ey Gl ae 


CONTINUE 


LF (LEV seO2 5) GO" te 
LP ( LEV Seo econ Re 
IF (LEVeEO: 3) G07 ro 
IF(IEV.EQ.4) GO TO 
EF (LEV@e ee 5) GCOs1o 


STOP 3004 
DOY 2 1 Ba iliae 
T( 1 yay ee) 
T(I,4) =YYY(I, 4) 
T (see Ce) 
TN ieoe— er ( Dean) 
CONTINUE 
GO TO 30 
DOm22 -l-176 
al ale 3) =VYY (TI, thy 
T(L, 4)=YYY ( bae 
2(L, 5) Shee Cie) 
AL (AG 6) =YYY (I, 6) 
CONTINUE 
GO °1e 30 


DOr 23 i= eG 
TL, 3] ae ls) 
T(1;4) 2yew 22) 
ECL, 5) =e 12S) 
TCL -6)'= vv (Ieee) 

CONTINUE 

GO TO 30 

DO 2452 = ic 


Wis) =YvYyY (ised) 
Pi, 2) =Vyvi £2 
T(1,5) =YYyY (1,3) 
TCI 6 y= vy 6) 

CONTINUE 

GO TO) 30 

DO 25 I=1,6 
Tey (le 
lita ana (| Tez) 
TI, 5) =vyy (1,3) 
T(I,6)=YYY(I, 4) 

CONTINUE 

CONTINUE 


NORMALIZATION OF THE CRITICAL EIGENVECTOR 
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a ha | 


CPerry aoe 


AA 


CALL NORMAL (T) 


INVERT TRANSFORMATION MATRIX 


CALL DLUD (6,6, TSAVE, 6, TLUD, IVLUD) 
pee4 I=1,6 

Bee. 2VLUD{(I) .EQ.0) STOP 3003 
CONTINUE 
@2EL DILU(6,6,TLUD, IVLUD, SVLUD) 
BOs I=1,6 

DO 9 J=1,6 

fey (1, J) =TLUD(I,J) 

CONTINUE 
CONTINUE 


SeeckK iInv(T) *A*T 


mw LT=1 

fee IMULT.EO.1) CALL MULT (TINV, ASAVE, T, A2) 
Geet iMULT.EO.0) STOP 3007 

P1=A2 (1,1) 

Bees (2,2) 

Pea (3,3) 

O=A2 (4,4) 

R=A2 (5,5) 

S=A2 (6, 6) 

Pere (21,*)P1,P2,P,90,R,8 


DEFINITION OF Nij 


WWHNHONNNNNFRFPRRRFR 
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Or) 


N33=TINV (3,3) 
N43=TINV (4, 3) 
N53=TINV (5,3) 
N63=TINV(6, i 
N14=TINV (1, 4) 
N24=TINV (2,4) 
Ne a= TINV (3, 4) 
N44=TINV (4, 4) 
NS54=TINV (5, 4) 
N64=TINV(6, 4) 
N15=TINV(1,5) 
N25=TINV (2,5) 
N35 = TENG 5) 
N45=TINV (4,5) 
N55=TINV(5,5) 
NGOS=TiAnNviesG, 5) 
N16=TINV (1, 6) 
NZo=TENY (2,7 6} 
N36=TINV (3, 6) 
N46=Tinyi 6) 
NSG=PTiNvis 6) 
N66=TINV (6, 6) 


DEFINITION OF Mij 


Mibt= Tas, 1) 
) 
iE) 
) 

M12=T(1, 2} 
) 
) 
) 


M13= T(1, 3) 
M23=T(2,3) 
Me j=20(5., 3) 
ae = 4a) 
M53=T(5,3) 


M64=T(6,4 
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Mga ni 2,5) 
Mog? (3,5) 
Mm45=T (4,5) 
horo—-T (1,6) 
Pio=T(2,6) 
M46=T (4, 6) 
Peto=1(5, 6) 
Proo=T (6, 6) 
Peewee (70,*) N11,N12,N13 
WRITE(71,*) N14,N15,N16 
WeLTE(72,*) N21,N22,N23 
Viet (73,*) N24,N25,N26 
oreeie(/4,*) N3l,N32,N33 
mete (75,*) N34,N35,N36 
WRITE (76,*) N41,N42,N43 
WRITE (77,*) N44,N45,N46 
Meee (78,*) N5l,N52,N53 
WReETE(79,*) NS4,N55,N56 
fore (80,*) N61,N62,N63 
WRITE (81,*) N64,N65,N66 
a K1=1./8.* ( (ALPHA2**3)+ALPHAO) / (ALPHA2 ) 
3 K2=3 .*A3-.5* (ALPHA2**2) /ALPHAO 
a K3=1./((B**2)*(DO**2) ) * (ALPHA2+A) * ( (ALPHA0O/ALPHAZ2) + (A**2) ) 
4 K=K1* (K2+K3) 
c 
3 Meant *, wh,k,ccrit 
‘a 
fame, (3 *D0**2) 
a 


1232567890123 4567890123456789012345678901234567890123456789012345 
5789012 
eee 2 1 **3-BL* (CCRIT**3 *KPSI**3*M41**3+KR* *3 *MO1**3+ 
Ky Ae 3 Mo ee 
br eGhimm.* 2 °KPol**2*KR*M41**2*M51+ 
Deen 2ekPol**2*KY*M41**2*M61+ 
fee ~KPOt~KR**2*M5S1* *2*M41+3 *CCRIT*KPSI*KY**2*M61**2*M41+ 
Peele ny = 2 *Mol~*2*M51+3 *KR**2*KY *MSa* ™2*M61+ 
eee) KrPol~KR~KY*M41*M51*M61) 
ue 22 ~ 5 —Bi* (CCRIT**3*KPSI**3 *M42**3+KR**3*M52**3+ 
& KY ~* 3 ~M62 **3+ 
& ee Cres ~ Zale > Lae ~KR*M42**2*M52+ 
& peecenli**2*Keem *25khy*M42** 2* M62 + 
fees Pe  *KPST *KR* * 2 *M52* *2*M42+3 *CCRIT*KPSI*KY* *2 *M62**2*M42+ 
& Peete 2 *MO2*~*2*M52+3*KR* *2* KY*M52* 82 *M62+ 
& peo eRET”“KPSI*KR*KY *M42*M52*M62) 
fee > M21 **2*M22-BL* (3 *CCRIT**3*KPSI**3*M41**2*M42+ 
& Die MS UY *2*M52+ 
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R — RR — KI —|F} -K 


RRR — — — RK —- -K 


SKY. *3 *M6 i? Zaheer 
3*CCRIT* *2*KPSI**2*KR* (M41 **2*M52+2 *M41 “M42 “Vie 
3*CCRIT**2*KPSI**2*KY* (M41**2*M62+2*M4 “Ma? *) oe 
3*CCRIT*KPSI*KR**2* (M51 **2*M42--2 “Mal 4*Merie eee 
3*CCRIT*KPSI*KY**2* (Mo1l**2*M42--2 “fie Meza ce ee 
3*KR*KY**2* (M61**2*M52+2 "Mo EP MiGZ = Mibal) = 
3*KR**2*KY* (M51**2*M624+2 7M5 Mao) 
6*CCRIT*KPSI**KR* KY* (M41*M51*M62+ (M41*M52+M42*Ma SiGe 
L24=3 *A3 *M21*M22**2-BL* (3*CCRIT**3 *KPSI**3*M41 *M4 2 
Baca e Siw ls Iii teya vy eyac: 
So Kye > MG MiG ae 
3*CCRIT* *2* KPSI**2*KR* (M42 **2*M51+2*M41 *Ma24eeee 
3*CCRIT* *2*KPSI**2*KY* (M42**2*M61+2*M4 1 *Ma27 ieee 
3*CCRIT*KPSI*KR**2* (M52* *2*M41 2 MS eM S72 ae 
3*CCRIT*KPSI*KY* *2* (M62**2*M41+2*M61 *M62 “Maze 
3*KR*KY* *2* (M62**2*M51+2 *MGleMiew aS 
3* KR* *2*KY* (M52 **2*M61+2 *MSiee Meio ee 
6*CCRIT*KPSI*KR*KY* (M42 *M52 *M61+ (M41*M52+M42*M51) *M62) ) 
se VPS) Se is) 
SA lo 5) ul i 
BS 3 (= ee cee idea Ze Mil 2 
b34=(=iee2 3) a AM 
RIT=NiIZ * L210 
RIZ=N1Z*b23+NiseGes 
R13 =N12*L24+N13*L34 
RI4Z=NIZei22 hee 2 
BZ L=N2Z2 LNZ3 oe 
RZ 2=N22e NZ oS 
R23=N22*L244+N23*L34 
R24=N22*L224+N23*L32 


K= (3 *R114+R134+R22+3*R24) /8 
WRITE (il,2001) wn, cena 
CONTINUE 
STOP 
FORMAT (’ ENTER MIN,MAX, AND INCREMENTS OF 
WN (stepstodivide) ’) 


FORMAT (’ ENTER DELTASAT. ’ ) 
FORMAT (’ ENTER A3’) 

FORMAT (’ ENTER CCRIT’ ) 
PORMAT@( ob 15 25) 

END 


SUBROUTINE DSOMEG (IJK, WR, WI, OMEGA, CHECK) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(6) ,WI (6) 
CHECK=-1.0D+25 
BOs inmate 
Pe VR) sbteeCHECK) GO tor 
CHECK=WR (TI) 
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IJ=1 
SONTINUE 
OMEGA=DABS (WI (IJ) ) 
Pavia)! CT.0.D0) LJK=1J 


Pe ( rs) .LT.0.D0) IJK=IJ-1 
RETURN 


ce ce ee eee eee ee es ee es es es es es es es es es es es es es es es es cs cs cc cc cr rc ee ei 
cece cee eee eee ee ees ee ee ee ees es ee es es es es es es es es es es es es es es es es es ee ee ss es es es ee ee i 


SUBROUTINE NORMAL (T) 

ivebicCIT DOUBLE PRECISION (A-H,0-2) 
DIMENSION T(6,6) ,TNOR(6, 6) 

Meee —0 (1,1) **2+T(1,2)**2 

Gee OMS (CFAC) .LE.(1.D-10)) STOP 4001 
morn (1,1)=1.D0 


mem 2, 1)=(T(1,1)*T(2,1)4+T(2,2)*T(1, = pang 
meres, 1) =(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 
TNOR(4,1)=(T(1,1)*T(4,1)4+T(4,2)*0(1, 2))/CFAC 
meee 5,1)=(T(1,1)*T(5,1)4+7T(5,2) *7T(1,2)) /CFAC 
mero, 1)=(T(1,1)*T(6,1)+7T(6,2)*T(1,2))/CFAC 
mrer(l,2)=0.D0 
TNOR (2, ee Ore AT 1) -Tt2 1) a) /CFAC 
TNOR (3 Pera: 2)*P(1, a) "3, LPT (12) YERAC 
TNOR (4, Pye (T (4,247 (1, ee (a> 1) A oekel 2) /CFAC 
Mem (5,2)=(T(5,2)*T(1,1) pr ct aie ri peieipe 
men (6,2)=(T(6,2) +P (1, et Gorey 12) / CF AC 
DO 1 I=1,6 

me 2 J=1,2 

T(I,J)=TNOR(I,J) 

CONTINUE 
CONTINUE 
RETURN 
END 


_— o_o ee eee eee eee eee ee ee ee ee ee ee ee ee cc ec eee eee ee ee es es es es es es es es ce eee ee 
—<—<— os oe eee ee ee ee ee ee ee ee ee ee ec cc ce ec ec cc ee ee 


SUBROUTINE MULT (TINV,A,T,A2) 

Poi CilT DOUBLE PRECISION (A-H,0-Z) 

PRMENSTON TINV(6,6),A(6,6),T(6,6),A1(6,6) ,A2 (6, 6) 
wee J=1,6 
eee J=1,6 

fwil,J)=0.D0 

a2 (1, B20 b0 


A(I,K)*T(K,J)+A1(I,J) 
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Wie WN 


Onl: 


CONTINUE 
CONTINUE 
CONTINUE 

DO 6t— aac 
DO 7 d=ie 

DO Sak —56 

AZiei ee) = 
CONTINUE 


CONTINUE 
CONTINUE 
DO Jie iao 

WRITE (*, 1201) (Ati) = ilo 
CONTINUE 
DO Ze — io 

WRITE (*, [OM tt Ged) =o 
CONTINUE 
DOwl) f—i1is 

WRITE (*, LOD) GAZ Crue oy 


CONTINUE 
WRETE (= Ohya ele) 
RETURN 
FORMAT (4E15.5) 
END 
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TINV (1, K) *Al(K oA 


APPENDIX B 
CRITICAL VALUE OF C FOR [X, X] BASIS 


PROGRAM NCCRIT 


234567 


HOPF BIFURCATIONS 
NOMOTO’S FIRST ORDER MODEL 
CALCULATIONS FOR CCRITICAL 
MiPLICIT DOUBLE PRECISION (A-H,0-Z) 
Peeper PRECISION K1,K2,K,LPSI,LY,LR,1Z,L, 
Maso, Ny, NR,NVDOT, NRDOT, NDRS, NDRB, KPSI, KRPKY, K3 


PeMENSION AMAT(6,6),FV1(6),1V1(6),YYY(6, 6) 

DIMENSION WR(6),WI(6),TSAVE(6,6),TLUD(6,6) ,IVLUD(6) ,SVLUD (6) 
DIMENSION ASAVE(6,6),A1(6,6) ,A2(6,6) 

OPEN (11,FILE=’CVWN1.RES’, STATUS=‘'NEW’ ) 

Ween (12,FILE=’CVWN2.RES’ , STATUS=’NEW’ ) 

Mean (13,FILE=’CVWN3.RES’, STATUS=’ NEW’ ) 

WRITE (*,*) ’ENTER MIN,MAX, AND INCREMENTS IN CCRIT’ 
READ (*,*) CMIN,CMAX,IC 

WRITE (*,*) ‘’ENTER MIN,MAX, AND INCREMENTS IN WN’ 
READ (*,*)  WNMIN,WNMAX, INCR 

WRITE (*,*) ‘ENTER WNO’ 

READ (*,*) WNO 

WEIGHT=435.0 


IZ =45.0 

L = ee) 
RHO i 

G = ae 
XG =0.0104 


MASS =WEIGHT/G 

MASS =MASS/(0.5*RHO*L**3) 
IZ =e (0), 5 *RHO*L**5) 
XG =XG/L 

YRDOT =-0.00000 

YVDOT =-0.03430 

YR =+0.00000 

ay =-0.10700 

YDRS <=+0.01241 

YDRB <2=+0.01241 

NRDOT =-0.00047 

NVDOT =-0.00000 

NR =-0.00390 

NV =-0.00000 


4} 


NBBS §—=-0-337>YDRS 

NDREB@e—+0. 252 -YRRs 

DH = (EZ -NRDOT) * (MAaSS— 1 pelea 

& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

Alls ((LIZ=NRDOT)*YV - (MASS* XG yRpel 1-4) Bu 


Al2Z=( (iZeNRDOT)* (-MASS+. vr) 
& MASS*XG-YRDOT) * (-MASS*XG+NR) ) /DH 
A21=( (MASS-YVDOT) *NV- (MASS2XCG=NV DOTS. ay Dr 


( 
( 
( 
A22=( (MASS-YVDOT) * (-MASS*XG+NR) - 
( 
( 
( 
( 


& MASS*XG-NVDOT) * (-MASS+YR) ) /DH 

Bll=({(IZ-NRDOT) *YDRS- (MASSexe-VE vot) Tike) ea 
B12=( (1Z-NRDOT) *YDRB- (MASS*XG=VRbOl) EEE Bn 
B21=( (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 


BZ = 5722 

Cl=(Al11*A22-AZ21T*Al2) *(AZl BAe 

C2=(A11+A22) * (A21*B1-Al1*B2)+B2* (All“AZZ-AZi Ale 
C=] (AZ ee) ee 


A=C1/C2 
B=C3/C2 
S 
SS. | IDs) 
ILMAX=1500 
C 
DO 1 If=) Nez 
@ 


WN =WNMIN+ (WNMAX-WNMIN) * (II-1) / (INCR-1) 
G DELLE awl 

ALPHA0=WN* *3 

ALPHA1=2 .15*WN**2 

ALPHA2=1.75*WN 

KPSI=-ALPHA1/B 

KY =-ALPHAO0/B 

KR =-(ALPHA2+A) /B 

GAMA0 =WNO* *3 

GAMA1=2 .15*WNO**2 

GAMA2=1.75*WNO 

LY=A+GAMA2 

LPSI=A* LY+GAMA1 

LR=A*LPSI+GAMAO0 
C2345 G7 

DO 2 J=1,2¢ 

CCRIT=CMIN+ (J-1) * (CMAX-CMIN) / (IC-1) 


G A MATRIX 
AMAT (1,1) =0.0 
AMAT (1,2) =1.0 
AMA Ge a a0 
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ooo © 


Oo PP @Go@@ 


AMAT (2,3 20 

Peer (2,4)=B*CCRIT*KPSI 
AMAT (2,5) =B*KR 

me ( 2,6) =B*KY 

AMAT (3,1)=1. 


R 
Been tl KPSt 
*KR 
-LR+B* KY 


510 
0 


DOWOrFool FORO OoCOOGOoO OSs 


AMAT (5,5 


AMAT (6,1 


St Be S 
Sore 


CALL RG(6,6,AMAT,WR,WI,0,2Z2ZZ,IV1,FV1, IERR) 
CALL DSTABL (DEOS, WR, WI, FREQ) 
B=CCRIT 

mee .Gr.1) GO TO 10 
DEOSOO=DEOS 

UO00=U 

LL=0 

SemrO 2 

DEOSNN=DEOS 

UNN=U 

PR=DEOSNN* DEOSOO 

ieee R.GT.0.D0) GO TO 3 
LL=LL+1 

See tLL.GT.3) STOP 1000 

ie — 0 

UO=U00 
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UN=UNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
UL=U0 

UR=UN 
DEOSL=DEOSO 
DEOSR=DEOSN 
U=(UL AUR Zou 
CCRIT=U 


0) 
“Cer Pa kWer St 
sss 


OoiororoaoaqoqcoqoqoorhWwW WwW or @a@aeaaqa al @ 


AMAT (5,2) 
AMAT (5,3) =LR 

BMAT (5,4) —B“CCRIT-KPS. 
AMAL (oy) —5 KR 


AMAT (5,6) =-LR+B*KY 
AMAT (6,1) =0.0 

AME GG 2) —O0 0 

AMAT (6,3) =LY 

AMAT (6,4) =1.0 

AMET Gyre) —Oaae 

AMAT (6,6) =-LY 


CALL RG(6,6,AMAT,WR,WI,0,22Z, IV1,FV1, TERR) 
CALL DSTABL (DEOS,WR,WI, FREQ) 

U=CCRIT 

DEOSM=DEOS 
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NM OR NM 


i= U 

PRL=DEOSL* DEOSM 
PRR=DEOSR* DEOSM 

meer rRG.GrT.0.D0) GO TO 5 
wO-—UL 

UN=UM 

Meo SsO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (UL-UM) 

mei Dir .GT.EPS) GO TO 6 
U=UM 

GO TO 4 

IF (PRR.GT.0.D0) STOP 3200 
UO=UM 

UN=UR 

Pao SO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (UM-UR) 

IF (DIF.GT.EPS) GO TO 6 
U=UM 

LLL=10+LL 

men lT=-U 

fee (LLL,*) CCRIT, WN 
UOO=UNN 

DEOSOO=DEOSNN 

CONTINUE 

CONTINUE 


O01 FORMAT (215) 
END 


—_ oe cee ee eee ee ee ee ee ee ee ee ee eee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee eee ee eee ee ee ee eee eee ee eee eee eee eee eee oe 
— oe eee eee cee ae cee cee cee ae cee cee ame cee ape cee ae cee ame gee oe eee eee eee ee eee eee cee cee cee oe ee ee oe ee ee ee ee ee ee ee ee ee ee ee ee eo ee ee eee ee oe 


SUBROUTINE DSTABL (DEOS, WR, WI, OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(6),WI(6) 
DEOS=-1.0D+20 
mor l I=1,6 
meee tWR(i).LT.DEOS) GO TO 1 
DEOS=WR (T) 
ee = 1 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 


— oe ee ce ee ees es es eee ees es es es es es es es es es es es es es es es es es es es es es es es es eee ee eee ee ee ee 
——— oe ie ee eee ee ee es es es es es es es ce ee es es es es es es es es es es es es es es es es cs ce cw se ee ee ee ee ie 
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APPENDIX C 
HOPF BIFURCATION PROGRAM FOR [X, X] BASIS 


PROGRAM KKPSI 
moPE BIFURCATIONS 
NOMOTO’S FIRST ORDER MODEL 


CALCULATIONS FOR K AND CCRITICAL IF KPSI CHANGES W/ C 


234567 
PIPLICIT DOUBLE PRECISION (A-H,O-Z) 


C) pik a a ad al > tal 


BeeoLeE PRECISION K1,K2,K,LPSI,LY,LR,1IZ,L, 


& MASS, NV,NR, NVDOT, NRDOT, NDRS, NDRB, KPSI,KR, KY, K3, 
& ieee. boo eo 32 oo, Loe et, bS2,853,L54, 
& Pee M2 MIS Mia MS MolG M21 ,M22>M23 ,M24,M25,M26, 
& M31,M32,M33,M34,M35,M36,M41,M42,M43,M44,M45,M46, 
& M51,M52 ,M53 ,M54,M55,M56,M61,M62,M63,M64,M65,M66, 
& Pee Nis ed Nies Nieo Newee N22 Ne 3 ,NZ4,N25,NZ6, 
& het NS2 N3oag-7NG4 Neos Noo 4h N42,N43,N44,N45,N46, 
& NSE, N52, Neg Ned Nese t5o, Nol, N62, N63 ,N6497N65,N66 
¢ 
MeMeENSLON AMAT(6,6),T(6,6),TINV(6,6) ,FV1(6),1IV1(6) ,YYY(6,6) 
DIMENSION WR(6) ,W1I (6), TSAVE(6,6) , TLUD (6,6) , IVLUD(6) , SVLUD(6) 
DIMENSION ASAVE(6,6),A1(6,6),A2(6, 6) 
meen (10,F ITLE=’CVWN1.RES’ , STATUS=’OLD’ ) 
OPEN (11,FILE=’AKPSI .MAT’ , STATUS=’ NEW’ ) 
C OPEN (12, FILE=’IREAL.MAT’ , STATUS=’NEW’ ) 
¢ SPEN (13,FILE=’RVALS.MAT’ , STATUS=‘NEW’ ) 
C OPEN (14, FILE=’ KKKKK.MAT’ , STATUS=’ NEW’ } 
WmeeGHitT=435.0 
7, =45...0 
L = 
RHO =] .94 
G ee 
'G =O. 0 104 


iWASS  =WEIGHT/G 

Mos = =MASS/(0.5*RHO*L**3) 
IZ a 5) * RHO*L** 5) 
XG =XG/L 

meeOoT =-0.00000 

memorT =-0.03430 

pale =+0 .00000 

pa =-0.10700 

YDRS =+0.01241 

YDRB =+0.01241 
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© (Sl Ba ae 


™ 


200 


205 
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NRDOT =-0.00047 

NVDOT eee o0e Co 

NR =-V. O0sgG 

NV =-0.00000 

NDRS  §]=0Fe3 7 abe 

NDRB =+02=2837>7YDR5 

DH =(IZ-NRDOT) * (MASS-YVDOT) - 


(MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *livee en 

Al2=((IZ-NRDOT) * (-MASS+YR) - 
(MASS*XG-YRDOT) * (-MASS*XG+NR) ) /DH 
A21=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
A22=( (MASS-YVDOT) * (-MASS*XG+NR) - 

(MASS*XG-NVDOT) * (-MASS+YR) ) /DH 
Blil=( (IZ-NRDOT) *YDRS- (MASS*XG=YREGi Nek mu 
B1l2=((IZ-NRDOT) * YDRE-(MASS=x*6-—YRpe a ieRe i br 
B21=( (MASS-YVDOT) *NDRS- (MASS* XG-NVDOT) * YDRS) /DH 
B22=( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 
Blab ed2 
BZ SH e7 lees 2 


WR ERE (* , 1004) 

READ Cz IWN 

INCR=IWN 

WRITE (*,*) ‘ENTER OBSERVER WN’ 
READ (*,*) WNO 

WRITE (*,1007) 

READ (*7*) 43 

DO is Dsat 


WRITE ao Gy) 
READ 447—) DO 
Cl=(A11*A22-A2 eal 2m (AZT * BIR a ee 


C2=(A11+A22) * (A21*B1-A11*B2) +B2* (A11*A22-A21*A12) 


C3 =-(AZ Et Bl -Ali eZ eae 
A=C1/C2 
BC 37/ Cz 
A3=0.0 


IBXO) EID TEIN Es 


READ {107 *) .GECRET Wi 
jenqakioie Ma Kiel 

ALPHAO =WN* *3 
ALPHA1=2 .15*WN**2 
ALPHA2=1 .75*WN 
KPSI=-ALPHA1/B 

KY =-ALPHAO/B 
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Pee =—(ALPHAZ+A) /B 
GAMA0=WNO* * 3 
GAMA1=2 .15*WNO**2 
GAMA2=1 .75*WNO 


LY=A+GAMA2 


LPSI=A* LY+GAMA1 
LR=A* LPSI+GAMAO 
2234567 


CE ee i oe ee an ae ae ee ae cet ei a ae ae at 


AMAT ( 


Hon ub th db Ww Wb We th Ut Wt tt th We ti wt tM a th ne th th a th th th th tnt 


CONTINUE 


too OO Ore 


EP 


Orrooi WwWwrootrProaradncnddcdoorRrRWWWo 


= 15 


510) 
PeckiT*KPSI 
aR 


KY 


FH 


So OO) ec 21 OG oa Oe 
Tp 
FH 


mer 1 ~KPST 


a ik 


iin “KRY 


au 
0 


Ki. 


CALL RG(6,6,AMAT,WR,WI1,1,YYY,1IV1,FV1,IERR) 
CALL DSOMEG(IEV,WR,WI, OMEGA, CHECK) 
WRITE (*,*) IEV 
WRITE (12,*) (WR(IREAL) , IREAL=1, 6) 
OMEGA0=OMEGA 
DO 5 I=1,6 
T(I,1)=YYY(I,IEV) 
T(I,2) =-YYY (I, IEV+1) 
CONTINUE 
IF(IEV.EO.1) GO TO 13 
IF (LEV 22052) 360, Tom 
IF (LEVeHO. 3) eCOnloers 
LF (TEV SEO. 4) "GOu Te 16 
IF (lEV SEG" 5) Comune 7 
STOP 3004 
DO 21 I=1,6 
T(I,3)=YYY(I,3) 
T(I,4)=YYY(I,4) 
T (ie yyy (1, Ss) 
Tio) a1 yy (176) 
CONTINUE 
GOueToO: 30 
DO@p2 2) L=iee6 
T(L, 3) =YYY (I, i) 
T(I,4)=YYY(I,4) 
ML, SysyYvvet 15) 
T( 1, 6) =v (1; 6) 
CONTINUE 
GO TO 30 
DO 235 0=i-6 
TS aeeaay | le 


( 

T(I,4) =YYY(I,2) 

a > 

T(I,6)=YYY(I, 6) 
CONTINUE 
eo ero aie) 

DO 24 I=1,6 

oe 3) =YYY(I,1) 

4) =YYY(I,2) 

oe S)=YYY(I,3) 

T(TI, 6) =YYY (I, 6) 
CONTINUE 
Cor Lem) 


DO 25 I=1,6 
T (es emery (1 
T(E, 4) avVYYi 1,2) 
io) SY a ees) 
(6) = Yoav 4 } 


ORs 
30 


CONTINUE 
CONTINUE 
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rere, er) 


AAA WwW”w 


MOeMAULZATION OF THE CRITICAL EIGENVECTOR 
CALL NORMAL (T) 
INVERT TRANSFORMATION MATRIX 


mom. L-1,6 
mes J=1,6 
TINV (I,J) = 
TSAVE (I,J) 
CONTINUE 
CONTINUE 
@aEL DLUD(6,6,TSAVE,6,TLUD, IVLUD) 
Bend I=1,6 
mot cVGUD(T) .EO.0) STOP 3003 
CONTINUE 
feet DILU(6,6, TLUD, [LVLUD, SVLUD) 
poems [=-1,6 
Pens, 6 
mei (ry) =TLUD (I,J) 
CONTINUE 
CONTINUE 


Oe, © 
aed ) 


Saeck Inv(T) *A*T 


IMULT=1 

IF (IMULT.EQ.1) CALL MULT(TINV, ASAVE,T,A2) 
fe (iMULT.EO.0) STOP 

Pn. (3,3) 

Q=A2 (4, 4) 

mn2 (5,5) 

S=A2 (6,6) 

fewer (21,*) P,O,R,S 


DEFINITION OF Nij 


( 

( 
et =TINV (3,1) 
N41=TINV (4,1) 
fo) —TINV (5,1) 
N61=TINV (6,1) 
fee —TINV (1,2) 
N22=TINV (2,2) 


iez—TINV (3,2) 
N42=TINV (4,2) 
Mee —TINV (5, 2) 
Bee —TINV (6,2) 
NIS=TINV (1,3) 


ol 


Oe 


N23 =PRENVAZ 
N33=Lilhh 3 
N43=TIA4 
N53=TineS 3 
N63=TiMViG, Ss 
NI4Z=TER ee 4 
N2Z24=TINV tz) 4 
N3 42h 
N44= 09 (4, 4 
N54=Ti5 , 4 
N64=TINV (6,4 
NiS= Piel, 5 
NZ5= = Ta | Ss 5 


O2 


@ 
C 
C 
c 
€ 
ce 
e 
C 
G 
€ 


Ri=1./8.* ( (ALPHA2**3)+ALPHAO) / (ALPHA2) 

K2=3 .*A3-.5* (ALPHA2**2) /ALPHAO 

K3=1./((B**2) * (DO**2) ) * (ALPHA2+A) * ( (ALPHAO /ALPHAZ2 ) + (A**2) ) 
K=K1* (K2+K3) 


meine ~*~, wn,k,ccrit 


mes/ (3*D0**2) 


M9456/8901234567890123456789012345678901234567890123456789012345 


eee 7 3 - BL (CCRIT**3*KPSI**3 *M41**3+KR**3 *M51**3+ 
MY *=3*M61**3+ 
ie 1 2 *KPol**2*KR*M4ieee2*M5i+ 
ea C@rimme = 2*KPST**2*Ky*M4i**2*Mo1+ 
eer tT * KPSI*KR* *2*M51**2*M41+3*CCRIT*KPSI*KY**2 *M61**2*M41+ 
chee 2 Mio ~*2*M51+3*KR* *2*KY*M51* *2 *M61+ 
peck in=KPSI*KR*KY*M41*M51*M61) 
fee M22**3—-BL* (CCRIT**3 *KPSI**3 *M42 * *3+KR* *3*M52**3+ 
Ky 43 Mio 2* *3 + 
pieenmen a. 72 ~KPSI**2*KR* M42 **2*M52+ 
Pi@@mer==2*KPSI**2*KY*M42**2*M62+ 
See? “Kol “~KR* *2*M52**2*M42+3 *CCRIT*KPSI*KY* *2*M62**2*M42+ 
Seed Ky a2 27 M62 **2*M52+3*KR**2*KY*M52**2*M62+ 
peck haKPol*KR*KY*M42*M52*M62) 
fete A> M2) **2*M22-BL* (3 *CCRIT**3*KPSI**3*M41**2*M42+ 
metre) M51 **2*M52+ 
SO NYewor Mol ~2*M62+ 
pie eo eo Kho 2 KRe (M41 **2*M52+2*M41*M42*M51) + 
eer tt o2 “KPol**2*KY* (M41**2*M62+2*M41*M42*M61) + 
feeeer Pf hrsi*KR**2* (M51%*2*M42+2*M51*M52*M41) + 
PeeORI@KPSil*KY**2* (M61**2*M42+2*M61*M62*M41) + 
Been. Ky **2* (M61**2*M52+2*M61*M62*M51)+ 
Beet ~~ 2*KY* (M51 **2*M62+2*M51*M52 *M61) + 
6*CCRIT*KPSI**KR*KY* (M41*M51*M62+ (M41*M52+M42*M51) *M61) ) 
eee 2 | *M22**2-BL* (3 *CCRIT**3*KPSI**3 *M41*M42**2+ 
& Sie * 3 *M51*M52* *2+ 
& Be 7 So “MO1*M62**2+ 


RMR RS KI RK 


PPR RR 


RF — —R KE RK RK HK KS 


oo 


3*CCRIT* *2* KPSI**2*KR* (M42**2*M51+2 *M41 "M47 eee 
3*CCRIT**2* KPSI**2*KY* (M42**274MG 12 * Mia ee eee 
3*CCRIT*KPSI*KR**2* (M52 **2 “Mal 2 Mas Meee 
3*CCRIT*KPSI*KY**2* (M62 **2*M414+2*M61 “M62 ~ Mae 
3*KR*KY**2* (M62* *2*M51+2*M6 aio Zee 
3*KR**2*KY* (M52**2*M6152 "Motes Zaitewe) 
6*CCRIT*KPSI*KR*KY* (M42 *M52 *M61+ (M41*M52+M42*M51) *M62) ) 
jl Rs iL // fo) ) Bs 
L32=(] oye 2 = 43 
33 = (=e ee 
3 4=( =) Zoe Mee 
C23456789012345678901234567890123456789012345678901234567690124a8 
L51=-BL* (CCRIT**3*KPSI**3*M41**3+KR**3 *M51**3+KY * “ge 
3*CCRIT* *2*°KPSI* *2°KR Ma ones 
3*CCRIT* *2*KPSI * te KY “M4 ile 
3 *CCRIT*KPSI*KR* *2*M51**2*M41+3 *CCRIT*KPSI * KY ®* *2 “MG eee 
3*KR*KY**2*M61* *2*M51+3 *KR** 2 “KY “Monee 2 eee 
6*CCRIT*KPSI*KR* KY “Mai Ms i Mews 
L52=-BL* (CCRIT* *3 *KPSI**3 *M42**3+KR* *3 *M52**3+KY ~ 93 
3 *CCRIT**2* KPSI** 27 KR=MAZ Ss a2 ile. 
S*CCRIT**2*KPSI**2"KY* Ma 2 = ort 
3 *CCRIT* KPSI * KR* *2*M52** 2*M42+3 *CCRIT*KPSI* KY * *2 *MG 7 
3*KR*KY **2*M62* *2*M52+3*KR~3Z “KY M57” Sea 
6*CCRIT*KPSI“~KR* KY *M422M5Z ica, 
L53=-BL* (3 *CCRIT**3*KPSI**3*M41**2*M42+3*KR**3*MS 
37KY* * 3 *M6de S22 PG 2h 
3*CCRIT* *2*KPSI* *2*KR* (M41**2*M52+2* M41 *Ma2 ~Meaeee 
3*CCRIT**2* KPSI**2*KY* (M41**2*M62-42 “M4 Ma one 
3*CCRIT*KPSI*KR* *2* (M51* *24M4 24-225 1 4MS2 ee 
3*CCRIT*KPSI*KY* *2* (M61 **2*M42-2 M6 > Moz 
3 *KR*KY**2*(M61**2*M52+2 *MG6 IG a ee 
3*KR**2*KY* (M51**2*M624+27M5 1 “Maze 
6 *CCRIT* KPSI*KR*KY* (M41*M51*M62+ (M41 *M52+M42 *M51) *M61) ) 
L54=-BL* (3 *CCRIT* *3 *KPSI* *3 *M41 *M42 **2+3* KR* *3 *M5 0 * Mazen 
SKY * = 3 *Me ieee = 22+ 
3*CCRIT* *2*KPSI* *2*KR* (M42**2*M51 +2 *Mat2Mo7 eee 
3*CCRIT* *2*KPSI**2*KY* (M42**2*Mol +2 7M4 Pia ope 
3*CCRIT*KPSI*KR**2* (M52 ** 2 Maa? “Mel Misc ee 
3*CCRIT*KPSI*KY* *2* (M62'** 2" Md 76 GZ 
3*KR*KY* *2* (M62* *2*M51 +2 +MiGIeaiiGZ aie 
3*KR*t2 *KY* (M52**2*M60+27Ms i" MeZaier,) 
6*CCRIT*KPSI*KR*KY* (M42 *M52*M61+ (M41 *M52+M42*M51) Maa 
RII=ENIP27i2TNI3* bse oe 
RIZ=N1276234-N13 73S hie bes 
R13=N12 *L24+N13 *L34+N15*L54 
R14=N12*L22+N13 *L32+N15*L52 
R21=N22 *L21+N23 *L31+N25*L51 
R2Z2=N22*L25+N23*b344N257 bos 
R23 =N22*L24+N23 *L34+N25*L54 
RZ24=-N2Z2*L22+4N23" boZaNZoe boc 


RRR RR KF KH 


RR FR KH 


R27 RI RE LF RI KI LK RR RK 


RRR —R RP KR RK 


o4 


Piel, ~*) Rl1,R13 pR22,R24 

an ReveR IT S+R22+3*R24)/S 

ewer (11,2001) WN,K,CERIT 
imcONTINUE 


(‘ ENTER NUMBER OF DATA’ ) 
1006 FORMAT (’ ENTER DELTASAT. ” ) 
(‘ENTER A3’) 


ee ee ee es ee ee ee es es es es ee ee es es es es es es es es ee ce cc cc ee ee ee eee 
eee eee see ee es es es es es ee es es es es es es es es es es es es es es es es es es es es es es ee eee es ees ee es ee ie 


SUBROUTINE DSOMEG (IJK,WR,WI,OMEGA, CHECK) 
IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 
DIMENSION WR(6) ,WI(4) 
CHECK=-1.0D+25 
DO 1 I=1,6 
IF (WR(I) .LT.CHECK) GO TO 1 
CHECK=WR (I) 
IJ=I 
CONTINUE 
OMEGA=DABS (WI (IJ) ) 
Mmemeiwilild).GT.0.D0) IJK=IJ 
fei Wl(lo).LT.0.D0) IJK=IJ-1 
RETURN 


SUBROUTINE NORMAL (T) 

mae otCIT DOUBLE PRECISION (A-H,0O-Z) 
PeMENSION T(6,6) , TNOR(6, 6) 

Meee =T(1,1)**2+T(1,2) **2 

mee DABS (CFAC) abE.(1.D-10)) STOP 4001 
moR(1,1)=1.D0 


eee ee, 2) (2,1) +T(2,2) *T(1,2)) /CFAC 
Meee) —~( b( 1,1) *P(3,1)+T(3,2) *T(1,2))/CFAC 
meee, oi (Mal) *T(4,1)4+7T(4,2) *T(1le2))/GFAC 
feet, t)=(T(1,1)*7T(5,1)+T(5,2) *T(1,2))/CFAC 
Meo, 1) —(1T(1,1)*T(6,1)+T(6,2)*T(1,2)) /CFAC 


MOR (1,2) =0.D0 
Megs, 2) —-1 2,2) *T(1,1)-T(2,1)*T(1,2))/CFAC 
ee, 2) = 13,2) *T(1,1)-T(3,1)*T(1,2)) /CFAC 
Mee 4, 2) =(7(4,2) *T(1,1)-T(4,1) *T(1,2)) /CFAC 
ee >, 2) *T(1,1)-T(5,1)*T(1,2))/CFAC 
mee, 2)—=(1(6,2) *T(1,1)-T(6,1)*T(1,2))/CFAC 
meet I=-1,6 

me 2 J=1,2 

at, J)=TNOR(I,J) 

CONTINUE 
CONTINUE 
RETURN 


D0 


Wir Ui 


OV ~] © 


— cc cee ce ce cr tec cece cer cere ce cece ccc mcr cece ce ce ee ee eee ee ee eee 
eel lc cc ee er A 


SUBROUTINE MULT (TINV,A,T,A2) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION TINV (6,6) ,A(6,6),1T(6, 6G) = stGr Go) AZ Gare) 
DO 1a 
DO: 2 J =i6 
Al (i jd)=02 DG 
AZ (i dO bo 
CONTINUE 
CONTINUE 
DOys. 26 
DO 47 0= ao 
DO 5 Kaisa 
Al (i, )=A(1;,K)*42(K, 0) Altes) 
CONTINUE 
CONTINUE 
CONTINUE 
DO "6 f— ieaG 
DO” 7 3 =18 ac 
DO 6.5 k= iG 
AZ iw) = 
CONTINVE 
CONTINUE 
CONTINUE 
bO 11 aaa 
WRITE. (7101) (25 ego J=1 6) 
Ika CONTINUE 
DO 12. iS 
WRITE (=3ebO)) (eo eae 
12 CONTINUE 
DO LOG 
WRITE (*7 10D) (AZ) een) 


TINV (1, K) *Al (Ko) -AZ eee 


10 CONTINUE 
WRITE (*, 20d) VAZ (ay 
RETURN 
deo FORMAT (4E15.5) 
END 
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